function result = gauss_int2d_rec(function_f, a, b, c, d, N_gauss_int2d)

n = int64(sqrt(N_gauss_int2d));
[gauss_points, gauss_weights] = gauss_points_weights_rec(n);

V = [a;
    b;
    c;
    d];

sum = 0.0;
if mod(n,2) == 0
    for i = 1:n/2
        for j = 1:n/2
            xi = gauss_points(i);
            eta = gauss_points(j);
            X1 = map_to_local(V, -xi, -eta);
            X2 = map_to_local(V, xi, -eta);
            X3 = map_to_local(V, xi, eta);
            X4 = map_to_local(V, -xi, eta);
            det_J1 = compute_jacobi_det(V, -xi, -eta);
            det_J2 = compute_jacobi_det(V, xi, -eta);
            det_J3 = compute_jacobi_det(V, xi, eta);
            det_J4 = compute_jacobi_det(V, -xi, eta);
            f_val1 = function_f(X1(1), X1(2));
            f_val2 = function_f(X2(1), X2(2));
            f_val3 = function_f(X3(1), X3(2));
            f_val4 = function_f(X4(1), X4(2));
            sum = sum + gauss_weights(i)*gauss_weights(j)*(f_val1*abs(det_J1) + f_val2*abs(det_J2) + f_val3*abs(det_J3) + f_val4*abs(det_J4));
        end
    end
else
    for i = 1:(n-1)/2+1
        xi = gauss_points(1);
        eta = gauss_points(i);
        if i == 1
            X = map_to_local(V, xi, eta);
            det_J = compute_jacobi_det(V, xi, eta);
            f_val = function_f(X(1), X(2));
            sum = sum + gauss_weights(1)*gauss_weights(i)*(f_val*abs(det_J));
        else
            X1 = map_to_local(V, xi, eta);
            X2 = map_to_local(V, xi, -eta);
            X3 = map_to_local(V, eta, xi);
            X4 = map_to_local(V, -eta, xi);
            det_J1 = compute_jacobi_det(V, xi, eta);
            det_J2 = compute_jacobi_det(V, xi, -eta);
            det_J3 = compute_jacobi_det(V, eta, xi);
            det_J4 = compute_jacobi_det(V, -eta, xi);
            f_val1 = function_f(X1(1), X1(2));
            f_val2 = function_f(X2(1), X2(2));
            f_val3 = function_f(X3(1), X3(2));
            f_val4 = function_f(X4(1), X4(2));
            sum = sum + gauss_weights(1)*gauss_weights(i)*(f_val1*abs(det_J1)+ f_val2*abs(det_J2) + f_val3*abs(det_J3) + f_val4*abs(det_J4));
        end
    end
    for i = 1:(n-1)/2
        for j = 1:(n-1)/2
            xi = gauss_points(i+1);
            eta = gauss_points(j+1);
            X1 = map_to_local(V, -xi, -eta);
            X2 = map_to_local(V, xi, -eta);
            X3 = map_to_local(V, xi, eta);
            X4 = map_to_local(V, -xi, eta);
            det_J1 = compute_jacobi_det(V, -xi, -eta);
            det_J2 = compute_jacobi_det(V, xi, -eta);
            det_J3 = compute_jacobi_det(V, xi, eta);
            det_J4 = compute_jacobi_det(V, -xi, eta);
            f_val1 = function_f(X1(1), X1(2));
            f_val2 = function_f(X2(1), X2(2));
            f_val3 = function_f(X3(1), X3(2));
            f_val4 = function_f(X4(1), X4(2));
            sum = sum + gauss_weights(i+1)*gauss_weights(j+1)*(f_val1*abs(det_J1) + f_val2*abs(det_J2) + f_val3*abs(det_J3) + f_val4*abs(det_J4));
        end
    end
end



result = sum;

end

function X = map_to_local(V, xi, eta)

N = [0.25 * (1 - xi) * (1 - eta), 0.25 * (1 + xi) * (1 - eta), 0.25 * (1 + xi) * (1 + eta), 0.25 * (1 - xi) * (1 + eta)];
X = N*V;

end

function det_J = compute_jacobi_det(V, xi, eta)

dN_dxi_deta = [-0.25 * (1 - eta), 0.25 * (1 - eta), 0.25 * (1 + eta), -0.25 * (1 + eta);
    -0.25 * (1 - xi), -0.25 * (1 + xi), 0.25 * (1 + xi), 0.25 * (1 - xi)];
J = dN_dxi_deta*V;
det_J = det(J);

end